function pvplot(grid,ubase,udofs,up,outfile,test_name)
% pvplot(up,udofs,outfile,test_name)
%
% Write a (u,p) solution up in a format readable by paraview (.vtk).
% The pressure is interpolated on the velocity nodes and velocity and
% pressure are represented with respect to the same grid. For P1 and
% P2 elements, the corresponding Paraview elements are used, while
% highr order elements are represented on an equivalent P1 mesh, since
% the is no native support in Paraview for them.

 fid = fopen(outfile,'w');
 
 % Write the (velocity) grid
 line = '# vtk DataFile Version 2.0';
 fprintf(fid,'%s\n',line);
 line = test_name;
 fprintf(fid,'%s\n',line);
 line = 'ASCII';
 fprintf(fid,'%s\n',line);
 line = 'DATASET UNSTRUCTURED_GRID';
 fprintf(fid,'%s\n\n',line);

 line = ['POINTS ',num2str(udofs.ndofs),' float'];
 fprintf(fid,'%s\n',line);
 for i=1:udofs.ndofs
   fprintf(fid,'%f %f %f\n',udofs.x(1,i),udofs.x(2,i),udofs.x(3,i));
 end

 elem_dofs = size(udofs.dofs,1); % number of dofs per element
 switch ubase.k % we have to reorder the nodes as in paraview
  case 1
   idx = [1 2 3 4];
  case 2
   idx = [1 2 3 4 5 8 6 7 9 10];
  otherwise
   error('Order not supported');
 endswitch
 line = ['CELLS ',num2str(grid.ne),' ',num2str(grid.ne*(1+elem_dofs))];
 fprintf(fid,'\n%s\n',line);
 for i=1:grid.ne
   fprintf(fid,'%i',elem_dofs);
   fprintf(fid,' %i',udofs.dofs(idx,i)-1);
   fprintf(fid,'\n');
 end

 switch ubase.k
  case 1
   cell_type = 10;
  case 2
   cell_type = 24;
  otherwise
   error('Order not supported');
 endswitch
 line = ['CELL_TYPES ',num2str(grid.ne)];
 fprintf(fid,'\n%s\n',line);
 fprintf(fid,'%i ',cell_type*ones(1,grid.ne-1));
 fprintf(fid,'%i\n',cell_type); % write the last one followed by \n

 % Write the data
 line = ['POINT_DATA ',num2str(udofs.ndofs)];
 fprintf(fid,'\n%s\n',line);
 line = 'VECTORS velocity float';
 fprintf(fid,'%s\n',line);
 for i=1:udofs.ndofs
   ii = (i-1)*3;
   fprintf(fid,'%f %f %f\n',up(ii+1),up(ii+2),up(ii+3));
 end
  
 fclose(fid);
return
